/*
 *  MatrixHandle.h
 *  NSE_SOLVER
 *
 *  Created by Patrick Zulian on 1/16/11.
 *  
 *
 */

#ifndef NSE_SOLVER_MATRIX_HANDLE
#define NSE_SOLVER_MATRIX_HANDLE

#include "Base.h"
namespace nsesolver {
	
	template< class MatrixType >
	class MatrixHandle {
		MatrixType &m;
		Range xRange, yRange;
		
        
        typedef typename MatrixType::Scalar Scalar;
        
	public:
		MatrixHandle(MatrixType &m, const Range &xRange, const Range &yRange)
		: m(m), xRange(xRange), yRange(yRange)
		{}
		
        
		operator MatrixType() const
		{
			MatrixType result(xmax(), ymax());
			
			for (int i = 0; i < xmax(); i++) {
				for (int j = 0; j < ymax(); j++) {
					result(i,j) = (*this)(i, j);
				}
			}
			return result;
		}
        
        
		
		inline int xmax() const
		{
			return xRange.getExtent();
		}
		
		inline int ymax() const
		{
			return yRange.getExtent();
		}
		
		inline int size() const 
        {
            return xmax() * ymax();
        }
        
		inline 
		typename MatrixType::Scalar & operator()(int i, int j) const
		{	
			return m(xRange.getFrom() + i, yRange.getFrom() + j);
            
		}
		
		const MatrixHandle & operator = (const MatrixType &matrix) const 
		{
			
			for (int i = 0; i < matrix.xmax(); i++) {
				for (int j = 0; j < matrix.ymax(); j++) {
					(*this)(i,j) = matrix(i, j); 
				}
			}
			return *this;
		}
		
		
		const MatrixHandle & operator = (const MatrixHandle &other) const
		{
			
			NSE_SOLVER_ASSERT( xmax() == other.xmax() );
			NSE_SOLVER_ASSERT( ymax() == other.ymax() );
			
			
			NSE_SOLVER_DEBUG(
                             if (this == &other) { NSE_SOLVER_ASSERT( false ); }
                             if (&m == &(other.m)) { NSE_SOLVER_ASSERT( false ); }
                             );
            
			for (int i = 0; i < xmax(); i++) {
				for (int j = 0; j < ymax(); j++) {
					(*this)(i, j) = other(i, j);
				}
			}
			return *this;
		}
        
        
        
        
        const MatrixHandle & operator = (const typename MatrixType::Scalar value) const
		{
			for (int i = 0; i < xmax(); i++) {
				for (int j = 0; j < ymax(); j++) {
					(*this)(i, j) = value;
				}
			}
			return *this;
		}
        
		
		friend 
		const MatrixType operator *(const MatrixType &m, const MatrixHandle &handle)
		{
			MatrixType mh = handle;
			return m * mh;
		}
		
		friend 
		const MatrixType operator -(MatrixType &m, const MatrixHandle &handle)
		{
			
			MatrixType mh = handle;
			return m - mh;
		}
		
		friend 
		const MatrixType operator +(MatrixType &m, const MatrixHandle &handle)
		{
			MatrixType mh = handle;
			return m + mh;
		}
		
        
        friend 
		const MatrixType operator +(const MatrixHandle &handle, const typename MatrixType::Scalar addend )
		{
			MatrixType m = handle;
			return m += addend;
		}
        
        
        
        const MatrixHandle & operator +=(const typename MatrixType::Scalar addend ) const
		{
			for (int i = 0; i < xmax(); i++) {
                for (int j = 0; j < ymax(); j++) {
                    (*this)(i, j) += addend;
                }
            }
			return  *this;
		}
        
 
		const MatrixType tranposed()
		{
			MatrixType mh = *this;
			return mh.transposed();
		}
		
		//Output
		friend
		std::ostream & operator <<(std::ostream &os, const MatrixHandle &handle)
		{
			
            for (int j = 0; j < handle.ymax(); j++) {
                for (int i = 0; i < handle.xmax(); i++) {
					os << handle(i, j) << "\t";
				}
				
				os << std::endl;
			}
			
			
			return os;
		}
		
		
	};
	
}

#endif
